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ABSTRACT 



C\| ■ We formulate the problem of the late accretion phase of the evolution of an 



isothermal magnetic disk surrounding a forming star. The evolution is described 



. by the six-fluid MHD equations, accounting for the presence of neutrals, atomic & 

C — ■ 

. molecular ions, electrons, and neutral, positively and negatively charged grains. 

Only the electron fluid is assumed to be attached to the magnetic field, in order to 
investigate the effect of the detachment of the ions from the magnetic field lines 
that begins at densities as low as 10 8 cm -3 . The "central sink approximation" 



is used to circumvent the problem of describing the evolution inside the opaque 



central region for densities greater than 10 11 cm -3 . In this way, the structure and 
evolution of the isothermal disk surrounding the forming star can be studied at 
late times without having to implement the numerically costly radiative transfer 



required by the physics of the opaque core. The mass and magnetic flux accu- 
mulating in the forming star are calculated, as are their effects on the structure 
& evolution of the surrounding disk. 

The numerical method of solution first uses an adaptive grid and later, after 
a central region a few AUs in radius becomes opaque, switches to a stationary 
but nonuniform grid with a central sink cell. It also involves an implicit time 
integrator; an advective difference scheme that possesses the transportive prop- 
erty; a second-order difference approximation of forces inside a cell; an integral 
approximation of the gravitational and magnetic fields; and tensor artificial vis- 
cosity that permits an accurate investigation of the formation and evolution of 
shocks in the neutral fluid. 



Subject headings: accretion - IS dust - magnetic fields - MHD - star formation 
- shock waves 
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1. INTRODUCTION 

Interstellar molecular clouds are known to be the sites of star formation. Typical sizes 
range from 1 to 5 pc, masses from a few tens to 10 5 M Q , and mean densities from 10 2 to 10 3 
cm -3 . Denser fragments (or cores) are formed within individual clouds. Approximately 1% of 
the mass of a cloud is in the form of dust particles (grains). These clouds are cold, with tem- 
peratures w 10 K, while their spectral lines have Doppler-broadened linewidths that suggest 
supersonic (but subAlfenic) internal motions. Polarization maps (Hildebrand et al. 1999; 
Rao et al. 1998; Lai et al. 2002) reveal large-scale, ordered magnetic fields in the plane of 
the sky. OH Zeeman observations (Crutcher et al. 1987; Crutcher et al. 1994) find magnetic 
field strengths between 10 and 150 /xG. 

Observations also reveal a degree of ionization Xi ^ 10~ 4 in cloud envelopes but < 10~ 7 
in dense cores (Casern et al. 1998; Williams et al. 1998; Bergin et al. 1999). In the deep 
interior of molecular clouds, cosmic-ray ionization dominates other mechanisms, but in cloud 
envelopes UV ionization is also important. A fraction of the dust grains becomes charged 
(negatively and/or positively depending on the density), and the grains thereby play an 
important role in the fragmentation of molecular clouds and star formation. 

A survey of protostars in the Orion Nebula (Beckwith et al. 1990) found that almost half 
of the star forming cores are surrounded by dusty disks. In addition, the geometry of proto- 
stellar fragments is often disklike over a wide range of lengthscales (Lay et al. 1994; Sargent 
et al. 1988; Kaifu et al. 1984). Nevertheless molecular clouds and their cores rarely exhibit 
rotation significantly higher than that of the background medium, and when they do, their 
angular velocities (typical core rotational rates are < 3 x 10~ 14 sec -1 ) imply centrifugal 
forces much too small to impose a disklike geometry through rotational support perpendicu- 
lar to the axis of rotation (Saito et al. 1996). These angular velocities are much smaller than 
those expected from angular momentum conservation from the initial galactic rotation (ini- 
tial angular velocity f2 ~ 10~ 15 s -1 and density n ~ 1 cm -3 ) (Goldsmith & Arquilla 1985). 
If angular momentum were to be conserved, centrifugal forces would not allow even the for- 
mation of interstellar clouds (Mouschovias 1991, eq. [2]). This is one way of stating the " 
angular momentum problem" of star formation. 

The presence of magnetic fields in interstellar molecular clouds poses another serious 
problem for star formation. The magnetic flux of an interstellar blob of mass comparable to 
a stellar mass is typically between two and five orders of magnitude greater than observed 
fluxes of even strongly magnetic young stars (Basri, Marcy & Valenti 1992; Wade et al. 1998). 
This suggests that flux loss takes place at some stage during star formation. 

Despite the fact that the critical mass for collapse of a pressure-bounded, isothermal, 
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nonmagnetic, nonrotating spherical cloud is only 5.8 M at density n = 10 3 cm -3 and 
temperature T = 10 K (Bonnor 1956; Ebert 1955; Ebert 1957), molecular clouds are not 
collapsing. No large-scale velocities characteristic of collapse (or expansion) are observed. 
Taking into account the contribution of the magnetic field, Mouschovias (1976a,b) calculated 
equilibria of initially uniform, cold clouds embedded in a constant-pressure external medium 
and threaded by a frozen-in magnetic field. These are oblate structures relative to the 
field lines that acquire an hour-glass shape. Mouschovias and Spitzer (1976) used these 
equilibrium states to determine analytically the critical mass-to-flux ratio above which the 
magnetic field cannot support a cloud against its self-gravity. Under typical molecular cloud 
conditions (B = 30 /xG and n = 10 3 cm -3 ), their expression yields a critical mass for collapse 
pa 5 x 10 2 M Q , consistent with observations. That critical mass-to-flux ratio is 20% smaller 
than the value calculated later by Nakano & Nakamura (1978) for an infinitessimally thin, 
uniform disk of infinite radius. 

Early on, supersonic turbulence was proposed either as a means of supporting the 
clouds or explaining the observed linewidths. But this scenario has the unattractive qual- 
ity of high energy requirements to sustain the motion, because such motions decay very 
rapidly (Field 1970). Since the observed linewidths are supersonic but subAlfvenic, they 
were attributed to large-scale, non-linear, oscillations, which are indistinguishable from long- 
wavelength, long-lived hydrodynamic waves (Mouschovias 1975). Such waves, because of 
their slow dissipation, have low energy requirements (Arons & Max 1975). Analyses of ob- 
served linewidths show a power-law relation between the velocity dispersion and the size of 
the observed objects (Larson 1981; Leung, Kutner & Mead 1982; Myers 1983), which was 
thought to be the signature of Kolmogorov turbulence. But this relation was explained quan- 
titatively as a consequence of virialized linewidths of hydromagnetic waves in magnetically 
supported, self-gravitating clouds (Mouschovias 1987a; Mouschovias & Psaltis 1995). 

The condition that the mass-to-flux ratio must exceed a critical value for a cloud to 
collapse is necessary but not sufficient. Other factors, such as external pressure and internal 
turbulent pressure, may aid or prevent the collapse, respectively. In nature, the mass-to- 
flux ratio of a molecular cloud as a whole is not expected to be very different from the 
critical value. If the ratio were much greater, the entire cloud would collapse very quickly. 
Velocities characteristic of this kind of collapse are not observed and the overall efficiency of 
star formation is generally low (M star /M c i oud < 0.1) (Frerking et al. 1987; Zhou et al. 1990; 
Andre et al. 1993). If the mass-to-flux ratio of a cloud as a whole were much smaller than the 
critical value, then the star formation timescale would be too long and clouds would rarely, if 
ever form stars. Molecular clouds are typically observed to be slightly subcritical (Carr 1987; 
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Loren 1989; Crutcher et al. 1993; Heiles et al. 1993; Crutcher et al. 1996; Shu et al. 1999). 1 

In the interiors of magnetically supported, subcritical clouds, stars can still form when 
self-gravity drives neutral matter through magnetic field lines (and the plasma) to form frag- 
ments (or cores), whose mass-to-flux ratio eventually exceeds the critical value (Mouschovias 
1977, 1978, 1979b). This process is referred to as ambipolar diffusion. Ambipolar diffusion 
was first proposed by Mestel and Spitzer (1956) as a means by which an interstellar cloud 
as a whole would reduce its magnetic flux. However, Mouschovias (1978, 1979b) pointed out 
that the essence of ambipolar diffusion is a redistribution of mass in the central flux tubes of 
a molecular cloud, not a loss of magnetic flux by a cloud as a whole. The ambipolar-diffusion 
timescale was shown to be three orders of magnitude smaller in the interior of a cloud than 
in the outermost envelope, where the degree of ionization is much greater, leading to a bet- 
ter coupling between neutral particles and the magnetic field (Mouschovias 1979b; 1987a,b). 
This suggests a self-initiated fragmentation (or core formation) on the ambipolar-diffusion 
timescale, which is less than 2 x 10 6 yr at a density of 10 4 cm~ 3 (see Mouschovias 1987a, eq. 
[12b]; or Mouschovias 1996, eq. [1]): 



where tq and r ni are the free-fall and neutral-ion collision timescales, respectively. This 
expression is essentially independent of geometry (see Mouschovias 1987a). The notion 
aspoused by proponents of turbulence-driven star formation, namely, that "rapid star for- 
mation" is inconsistent with ambipolar-diffusion, self-initiated star formation, has no merit. 
The issue is discussed quantitatively in another publication (Tassis & Mouschovias 2004c). 
The small degree of ionization in cloud interiors also allows ambipolar diffusion to dump 
the hydromagnetic waves or magnetic turbulence (thus rendering them irrelevant to star 
formation) over a characteristic lengthscale comparable to the critical thermal lengthscale 
(~ Jean's length), thereby also leading to the thermalization of linewidths in agreement with 
observations by Myers & Benson (1983) and Myers, Linke & Benson (1983). This point of 
view is substantiated by multifluid MHD numerical calculations in rectilinear (Paleologou 
& Mouschovias 1983; Mouschovias, Paleologou & Fiedler 1985), cylindrical (Mouschovias 



1 Crutcher (1999), after compiling measurements from a number of molecular clouds, concluded that the 
typical cloud mass-to-flux ratio is a factor of two grater than the critical value. However, both he and Shu et 
al. (1999) recognize that the actual magnetic field strength must on average be a factor of two greater than 
the measured value, since only the line of sight component is measured by the Zecman effect. In addition, the 
average column density measured is by a factor of two greater than the actual relevant value (perpendicular 
to a thin disk) due to geometric effects. Applying these corrections, Shu et al. (1999) concluded that the 
clouds are typically subcritical. 
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& Morton 1991, 1992a,b) and axisymetric geometries (Fiedler & Mouschovias 1992, 1993; 
Ciolek & Mouschovias 1993; Morton, Mouschovias & Ciolek 1994; Ciolek & Mouschovias 
1994; Desch & Mouschovias 2001; Eng & Mouschovias 2004). 

For a cloud as a whole and for a core during its ambipolar-diffusion-controled, quastistatic 
(or subcritical) phase of contraction, magnetic braking is very effective and keeps the core and 
the cloud essentially corotating with the background (Mouschovias 1977, 1979a; Mouschovias 
& Paleologou 1979, 1980; Mouschovias & Morton 1985a,b; Basu & Mouschovias 1994, 
1995a,b). Hence the centrifugal forces have a negligible effect on the evolution of the con- 
tracting core. Magnetic braking is responsible for allowing clouds to form in the first place 
and for reducing the angular momenta of cloud cores to their observed low values, as well as 
for resolving the angular momentum problem of star formation during the early, isothermal 
stage of contraction, as required by observations (see review by Mouschovias 1981, §6) 

After the central mass-to-flux ratio exceeds the critical value, the now thermally and 
magnetically supercritical core begins to contract dynamically, while the envelope remains 
magnetically supported. The first axisymmetric calculations to follow both the quasistatic 
and the dynamic phase of core formation and contraction (Fiedler & Mouschovias 1992, 1993) 
found that force balance along field lines is established rapidly and maintained even while 
dynamical contraction takes place perpendicular to the field lines. A disklike geometry is 
established early and maintained throughout the evolution of the protostellar fragment. This 
result was exploited by Ciolek & Mouschovias (1993, 1994) and Basu & Mouschovias (1994, 
1995a,b), who assumed force balance along field lines and integrated the equations analyti- 
cally over the vertical structure of the disk, thus reducing the dimensionality of the problem 
for computational purposes without losing any of the essential physics of the problem. 

These axisymmetric simulations have shown that, even during the dynamical phase of 
contraction, the infall of a core does not evolve into free fall (see review by Mouschovias 
1995). Even at high central densities (n = 10 10 cm -3 ), the maximum infall acceleration 
does not exceed 30% that of gravity, and the magnetic force dominates the thermal-pressure 
force almost everywhere in the supercritical core except in a small central region. During 
the dynamical contraction phase, both magnetic flux and angular momentum tend to get 
trapped inside the core, with a consequent increase of the magnetic field and angular velocity 
with gas density as B c oc p l J 2 and Q c oc p l J 2 (Mouschovias 1976b, 1979a, 1989; Fiedler & 
Mouschovias 1993; Basu & Mouschovias 1994). 

Millimeter and sub-millimeter continuum observations (Ward-Thompson et al. 1994; 
Andre et al. 1996; Bacmann et al. 2000) yield density profiles for starless cores that are in 
agreement with those predicted by the ambipolar diffusion models. Specific models con- 
structed for the Barnard 1 cloud (Crutcher et al. 1994) and L1544 (Ciolek & Basu 2000) 
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found the core properties predicted by the ambipolar diffusion model to be in excellent 
agreement with the observed values. In addition, Benson, Caselli, & Myers (1998) measured 
ion-neutral drift speeds consistent with the theoretical predictions. 

The leftover magnetic flux at the end of the early, isothermal stages of contraction 
(n < 10 10 cm~ 3 ) exceeds observed stellar magnetic fluxes and must therefore be either 
dissipated or somehow excluded from the matter destined to give birth to stars. Approximate 
theoretical estimates argued that Ohmic dissipation is effective in destroying the magnetic 
flux at densities above 10 10 cm -3 (Pneuman & Mitchell 1965; Nishi, Nakano & Umebayashi 
1991). Recent dynamical calculations by Desch and Mouschovias (2001), which follow the 
evolution up to densities n = 2 x 10 12 cm -3 , show that ambipolar diffusion is much more 
effective than Ohmic dissipation in decoupling the magnetic field from the neutral matter. 
(Earlier work overestimated the e — H 2 elastic collision cross section, and hence the ohmic 
resistivity, by a factor of 100; see review by Mouschovias 1996). 

The evolution timescale of the dynamically contracting supercritical core at the end of 
the above simulations is less than 100 yr; a protostar in hydrostatic equilibrium is expected 
to form soon thereafter. There have been several studies of the accretion phase following the 
formation of a protostellar core. Tomisaka (1996) performed a two-dimensional, axisymet- 
ric, isothermal simulation and used the "sink cell method" to study the accretion phase. He 
started from a supercritical initial condition and assumed flux-freezing throughout the sim- 
ulation. Shu and collaborators extended the singular isothermal sphere similarity solution 
(Shu 1977). Galli & Shu (1993a,b) included a spatially uniform, dynamically weak magnetic 
field, while Li & Shu (1997) studied the accretion of a singular isothermal disk with a uniform 
mass-to-flux ratio, assuming the flux to be frozen in the neutral fluid. In addition, in these 
studies the authors always assume that the accretion onto the protostellar core starts from a 
static singular configuration. This assumption is inconsistent both with the aforementioned 
analytical and numerical studies of the evolution before the formation of a hydrostatic proto- 
stellar core, which find supercritical cores collapsing dynamically, and with observations that 
find significant infall velocities in prestellar cores (Tafalla et al. 1998; Williams et al. 1999). 
Safier, McKee & Stahler (1997) used a spherical model to study the collapse of both a 
prestellar core as well as the accretion phase after the formation of a protostellar core. They 
neglected the thermal pressure and the magnetic tension force and they included ambipolar 
diffusion in a phenomenological way without solving the induction equation for the self- 
consistent evolution of the magnetic field. Their spherical model was extended by Li (1998a, 
1998b) to include the thermal pressure and the induction equation. However, these models 
do not satisfy the solenoidal condition on the magnetic field and cannot account for the 
magnetic tension force. 
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Taking into account the decoupling of the magnetic field from the matter Li & McKee 
(1996) suggested that Ohmic dissipation would halt the accretion of flux onto the protostar 
and give rise to a hydromagnetic shock. However, Ciolek & Konigl (1998), extending the 
previous numerical studies of Mouschovias and collaborators to study the later accretion 
phase of star formation by using the "central sink cell" approach and assuming the magnetic 
flux to be frozen in the ion fluid, found that ambipolar diffusion is sufficiently effective at 
much lower densities (or much larger radii) to halt flux advection and drive a hydromag- 
netic shock, independent of the effect of Ohmic dissipation. This hydromagnetic shock is 
also present in the self-similar solution obtained by Contopoulos, Ciolek & Konigl (1998). 
Starting from a dynamical initial state (nonzero infall velocity) and assuming force balance 
in the z-direction as well as a power-law scaling of the density of ions with respect to the 
density of neutrals, they were able to obtain a similarity solution which is very close to the 
numerical results of Ciolek & Konigl (1998). 

In this paper we formulate the problem of the late accretion phase of an isothermal 
magnetic disk surrounding a protostar by taking into account important physics ignored by 
previous calculations: (1) the decoupling of the ions from the magnetic field lines, which 
occurs at densities above 10 8 cm -3 (Desch & Mouschovias 2001) and which leaves the mag- 
netic field frozen in the much less massive and much more tenuous electron fluid alone; (2) 
the chemical and dynamical effects of the positively-charged grains, the abundance of which 
becomes significant in the same density regime. The results of this calculation are discussed 
in a companion paper (Tassis & Mouschovias 2004; hereafter Paper II). 

In §2 we describe the model cloud. In §3 we discuss the six-fluid MHD equations 
appropriate for the physical system under study, and in §4 the thin-disk approximation 
employed in our numerical code. In §5 we state our initial and boundary conditions, and in §6 
we discuss in detail the implementation of the central sink approximation. The dimensionless 
equations solved numerically and the method of solution are described in §7. We conclude 
in §8 by discussing the most important differences of our formulation from previous work 
applied to similar problems. 

2. MODEL CLOUD 

We consider a self-gravitating, magnetic, weakly-ionized model molecular cloud con- 
sisting of neutral particles (H2 with 20% He by number), ions (both molecular HCO + 
and atomic Na + , Mg + ), electrons, neutral grains, singly negatively- charged and positively- 
charged grains. Following Desch and Mouschovias (2001), the abundances of charged species 
is determined from the chemical reaction network shown in Table 1. (For a detailed dis- 
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cussion of the chemistry, see Ciolek & Mouschovias 1993, 1995). The positively-charged 
grains cannot be neglected at densities grater than 10 9 cm -3 , since their abundance becomes 
comparable to that of negative grains (see Fig. 2 in paper II). For simplicity, we consider 
spherical grains of uniform size. Desch and Mouschovias (2001) used a chemical model that 
accounted for grains of different sizes and found that the effect of the grain size distribution 
on the evolution of the cloud is minimal. In the case of collisions of ions (molecular or 
atomic) with grains, we assume that the ions do not get attached to the grains, but that 
they get neutralized, with the resulting neutral particle escaping into the gas phase. Thus 
the total abundance of metals as well as the HCO abundance remain constant. Grain growth 
and depletion of elements by attachment onto grains are not considered in this paper. 

The axisymmetric model cloud is initially in an exact equilibrium state with self-gravity 
and external pressure being balanced by internal magnetic and thermal-pressure forces. We 
follow the ambipolar-diffusion-initiated evolution up to a central density of 10 11 cm" 3 , below 
which isothermality is a good approximation. To study the later stages of the evolution of the 
disk surrounding the forming star, we employ the "central sink cell" method, originally used 
by Boss & Black (1982) and more recently employed by Ciolek & Konigl (1998). The radius 
of the central sink defines the inner boundary r = -Ri nner of the computational region. The 
use of a central sink makes the study of the dynamics of the collapsing protostellar fragment 
possible, even after the formation of an opaque core, by fixing the value of i?i nne r at a radius 
beyond which the density is low enough (n n < 10 11 cm -3 ) that isothermality remains an 
excellent approximation. Computationally taxing radiative transfer becomes unnecessary. 

3. SIX-FLUID MHD DESCRIPTION OF MAGNETIC STAR FORMATION 

We extend the multifluid formalism described by Ciolek and Mouschovias (1993) to 
include the positively charged grains and to relax the assumption of flux-freezing in the ions. 
Only the electrons are assumed to be attached to the magnetic field lines. The equations 
governing the behavior of the six-fluid system (neutrals, electrons, ions, negative, positive, 
and neutral grains) are: 



dp, 



+ V • ( Pn v n ) = 0, 



(la) 



dt 




(lb) 




(lc) 
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= -en e (E + -j x B) + F en , (Id) 

= em(E + — x B) + Fin, (le) 

c 

= -en g _ (EAxB)l F g _ n + F g _ g0)inel , (If) 



= en g+ (F + x B) + F g+n + F g+g0 , inel , (lg) 
= F gon + F gog) i ne i + F gog+! ; ne i, (lh) 

VxB = — j, (li) 

c 

j = ein-^i + n g+ v g+ - n c v c - n g _v g _), (lj) 

— = -c(V x E), (Ik) 
V • flf = -4ttG Pq , (11) 

Pn = PnCl (lm) 

The quantities p s , n s and v s denote the mass density, number density, and velocity of species 
s. The subscripts n, i, e, g_, g + and g refer to the neutrals, ions, electrons, negatively- 
charged grains, positively-charged grains, and neutral grains, respectively. The symbols g, 
E and B denote the gravitational, electric and magnetic fields, respectively. The magnetic 
field satisfies the condition V • B — everywhere at all times. 

The frictional force (per unit volume) on species s due to elastic collisions with neutrals 
is given by 

F sn = -K - v s ), s = i, e, g_, g + , g , (2) 

where the mean collision time, accounting for both s - H 2 and s - He collisions, is 

T S H 2 1 ^H 2 + m s 
Tsn= = j r , S = l,e,g_,g+,g , 

and the quantity 



(3) 



a sRc = 1.14 for s = i, (4) 
= 1.16 for s = e, (5) 
= 1.28 for s = g_,g+ or g (6) 

is the factor by which the presence of He reduces the slowing-down time of species s with 
respect to s - H 2 collisions alone (see Mouschovias 1996, §2.1). The mean collisional rate 
between ions of mass and hydrogen molecules of mass mu 2 is (crw)iH 2 = 1.69 x 10~ 9 cm 3 s _1 
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for HCO + — H 2 collisions, and almost identical to this value for Na + — H 2 and Mg + — H 2 
collisions. The value of the mean collisional rate between electrons and H 2 is {aw) e n 2 = 
1.3 x 10~ 9 cm 3 s _1 . In calculating these collisional rates, the Langevin approximation is used 
for ion-neutral collisions but not for electron-neutral collisions, for which the electron spin 
is important (Mouschovias 1996). The grain-H 2 collisional rate is given by 

(aw) gU2 = 7ra 2 (8A;BT/7rmH 2 ) 1/2 . (7) 

This expression is valid only if the velocity difference between a grain and a hydrogen 
molecule is smaller than the sound speed in the neutrals (Ciolek & Mouschovias 1993). Oth- 
erwise, one has to use 

((Tw)gH 2 = ira 2 \v n — v g \. (8) 

The quantity F 7 s,me\(k) i n equations (If) - (lh) is the force (per unit volume) on grain 
fluid 7 due to the conversion of dust particles of fluid 5 into dust particles of fluid 7, with 
mass gain 8m 1 , (note that m 7 = ms) via an inelastic process k, where k can be any one of 
the processes listed in Table 1 that involve grains. The momentum transferred to fluid 7 
from fluid 5 is 

Ap = \Sm y \(v s - u 7 ); (9) 

and we can write 



■^* 1 7i5,inel(fe) m 7 



<9n 7 



()1 (vs-v 7 ), (10) 

where the time rate of change of the number density of 7 due to the inelastic process (k) is 
given by 



"7 



dt 



= n v n 5 ((Tw) {k) 
(*) 

(id 



7~<5?7,mel 



and r^inei is the timescale for a particle 5 to find a particle 77 and be converted into 7. Thus, 
in general, since m 7 = ms 

Ps 

Fy8,incl(k) = (VS - V 7 ). (12) 

7~<5?7, inel 

Furthermore, according to Newton's third law, the fluid S will also experience an equal and 
opposite force 

-FYy,incl(fc) = —Fy5,mcl(k) = (t> 7 — V$). (13) 

^^incl 

The relevant timescales due to these inelastic processes are r g o e ,inei = (^e^ego) -1 ; T g oi,inei = 

(TiittjgCl) 5 7g_i,inel (^i^ig-) j 7~g + c,incl (j^e^eg+) 5 ^g_g+,inel (^g+ ®g-g+ ) ' ^g+g-^nel 
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(n g _ai g _) 1 , where the reaction rates a cg o, «i g o, «cg + and a g _ g+ are given in Appendix 
A. 

In the force equations of the electrons, ions, and grains, the acceleration terms have been 
neglected because of the small inertia of these species. Mouschovias, Paleologou & Fiedler 
(1985) included the acceleration term for the plasma and showed that the plasma reaches 
a terminal drift velocity very fast. Likewise, the thermal-pressure and gravitational forces 
(per unit volume) have been dropped from the force equations of all species other than the 
neutrals because they are negligible compared to the electromagnetic and collisional forces. 
The inelastic momentum losses by the electron and ion fluids due to attachment onto grains 
and neutralization are negligible compared to the momentum loss due to elastic collisions, 
and they have been omitted from the force equations (Id) and (le). 

Finally, since we constrain our study to the region characterized by densities smaller 
than 10 11 cm~ 3 , isothermality is a good approximation, and C n = (k^T / /jmn) 1 ' 2 is the 
isothermal speed of sound in the neutrals, &b the Boltzmann constant, and \i the mean mass 
per neutral particle in units of the atomic- hydrogen mass, mn- For the early stages of star 
formation or for the isothermally infalling gas outside the optically thick inner region, the 
temperatures are much smaller than the grain sublimation temperature ~ 1500 K, and a 
mass continuity equation (lb) for the grains (charged + neutral) is appropriate. 

Altogether, then, we have a system of 13 equations, (la)-(lm), which contain 17 un- 
knowns (p n , P, E, B, j, g, v 

m ) ^g+) ^go> Pc? Pi? Pg- ? Pg+ 5 Pgo)' close the system, 

the densities of electrons, ions and charged grains (n e , n i? n g „, and n g+ ) are calculated from 
the equilibrium chemical model described in Appendix A. 

Solving the electron force equation for the electric field, we obtain the generalized Ohm's 

law 

E = -^xB + ^ Vc - V «\ (14) 

The second term on the right-hand side can be neglected for densities < 10 14 cm -3 , reducing 
equation (14) to the simple form 

E = -—xB. (15) 

c 

This is inserted in Faraday's law of induction (eq. [Ik]) to obtain the flux- freezing equation 
in the electrons 

OB 

— = Vx(v c xB). (16) 
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4. GOVERNING EQUATIONS IN THE THIN-DISK APPROXIMATION 

Numerical simulations of the ambipolar-diffusion induced evolution of axisymmetric, 
isothermal, molecular clouds have shown that force balance is rapidly established along field 
lines and is maintained throughout the evolution of the model cloud, even during the dynamic 
phase of contraction that follows the formation of a thermally and magnetically supercritical 
core (Fiedler & Mouschovias 1992,1993; Desch & Mouschovias 2001). Following the approach 
of Ciolek & Mouschovias (1993), we take advantage of these results and model the cloud as 
an axisymmetric thin (but not infinitesimally thin) disk, with its axis of symmetry aligned 
with the ^-axis of a cylindrical polar coordinate system (r, 0, z) and with force balance along 
the field lines maintained at all times. The sense in which the thin-disk approximation is 
used is that the radial variation of any physical quantity is small over a distance comparable 
to the local disk half-thickness. The upper and lower surfaces of the disk are at z = +Z(r, t) 
and z = —Z(r, t), while the radius R of the disk satisfies the condition R^> Z. The magnetic 
field has the form 

B(r, z, t) — zB Zieq (r, t) for \z\ < Z(r,t), (17a) 

B(r,z,t) = zB z (r,z,t) +fB r (r,z,t) for \z\>Z(r,t). (17b) 

Inside the disk the magnetic field has only a z-component, equal to the field in the equatorial 
plane B z ^ cq . The r-component of the field is antisymmetric about the equatorial plane, i.e., 
B r (r, z) = —B r (r, —z). Far away from the cloud the magnetic field is assumed to be uniform, 
perpendicular to the plane of the disk and constant in time. The normal component of the 
magnetic field is continuous across the upper and lower surfaces of the disk. The unit vector 
normal to the upper and lower surfaces has both z and r components (i.e., the surface is not 
parallel to the equatorial plane). The model cloud is embedded in an external medium of 
constant thermal pressure P cxt . 

For computational simplicity, we reduce the mathematical dimentionality of the prob- 
lem by integrating the equations analytically over z, assuming ^-independence of all physical 
quantities. This is an excellent approximation since thermal-pressure forces smooth out den- 
sity gradients over lengthscales smaller than the critical thermal lengthscale (see Mouschovias 
1991), which is comparable to Z(r,t). The thin-disk equations governing the model cloud 
are 

r+Z(r,t) 



p D dz = 2p n (r,t)Z(r,t), (18a) 

■Z(r,t) 

= 0, (18b) 



-Z(r,t) 

da n 1 d{ra n v ^ r ) 
dt r dr 



djXgVn) | 1 d[ra Q {x g ~v g ~, r + Xg+v g +,r + X g o%),r)] _ Q , lgc v 
dt r dr ' 
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Pn (r,t)C 2 = P cxt + ^Gal(r,t), 



d(a n v n:r ) 1 d{ra n vl r ) 



dt 



r dr 



n 2 d °n 
— — O c 



'eff' 



dr 



+ a n g r + F l 



mag,n 



vr 2 [3P ext + (7r/2)G^] ^ 2 



2 n [Pext + (vr/2)G^] 



(18d) 
(18e) 
(18f) 



mag,r 



18Z 

2~dr~ 



dZ 



B r ,z + 2-B 2 , eq ( -B r ^— J + ( B T)Z 



dr 



M(r,r') = ^- 
dr 



POO 

g r (r) = 2nG dr'r'a n (r')M(r,r'), 
Jo 

POO 

/ dkJ (kr)J (kr') 
Jo 

POO 

B r ,z(r) = - / dr'r'[B zm {r') - B rci ]M(r, r% 
Jo 

dB ZtCq 1 d(rB Zfiq v c , r ) = Q 
dt r dr 



(18h) 



2 d 


-><(-)] 




7r dr 


r> V r >/. 



We = f n + 



A mag ' ne 



<7 n A e 



-V n + 



6i + 1 0; + 1 



i e ; „ 

-V n + 



(18n,o) 



(18p,q) 



u g- = 



7 g+ 



1 O 



i 



(18i) 

(18j) 
(18k) 

(181,m) 



e _ + l e _ + l 



-V n + 



e+ + i ©g+ + i 



V e , v g-,(f> "g-,? 



^e, Wg+,0 — — g+,, 



0„ 



g -,^ + 1 e g _,^ + i 
e 



V n + - S+ ^^ e,. 







-V n 



gO 



0gO + 1 ©HO + 1 



go 



^e, %),</> — "gO,c 



1 



©g+,0 + 1 0g+,</- + 1 



-V n 



-V Q + 



0gO,^ + 1 0gO,4> + 1 

0e .<h 



18r,s) 

e^ + r u ' e^ + i'V (18t) 

All symbols introduced in equations (18k)-(18t) are defined in Appendix B. A detailed deriva- 
tion of the first nine thin-disk equations, (18a)-(18i), is given in Ciolek & Mouschovias (1993). 
Introducing the effective sound speed, C e s , we account for the contribution of the radial force 
exerted by the external pressure P cxt on the upper and lower surfaces of the disk, since the 
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latter are not horizontal. In the expression for the total radial magnetic force (per unit area), 
-Pmag,r, includes the magnetic tension force (the B r terms) acting on the surfaces of the disk. 
The r-component of the gravitational field is calculated from Poisson's equation for a thin 
disk (18h). J is the zeroth-order Bessel function of the first kind, K is the complete elliptic 
integral of the first kind (Abramowitz & Stegun 1965), r< = min[r, r'], and r> = max[r, r']. 

Equations (181-t) give the velocities of the ions and the three grain species in the r and 
directions, and the 0-component of the electron velocity in terms of the radial velocities of 
the neutrals (v n ) and electrons (i> e ; eq. [18k]). They come from the solution of the steady- 
state force equations of the species and Ampere's law (eqs. [ld]-[lh] and [li], where j in 
eq. [li] is given by eq. [lj]). These are algebraic equations, but they have coupled r and 
(j) components due to the presence of the Lorentz terms (no z-components since we assume 
^-independence for all quantities inside the model cloud). All the radial velocities, except 
that of the electrons, have the form 

Bs -v e + ^—rV n (19) 



3 e s + i c e s + i 11 

(see Appendix A), where Q s is the indirect attachment parameter of species s : for @ s ^> 1 
the r-velocity of species s is approximately equal to that of the electrons, and the species 
is well attached to the field lines, while for Q s <C 1 the species is detached and is falling in 
with the neutrals (v s ~ v n ). 



5. INITIAL AND BOUNDARY CONDITIONS 

The model cloud is initially in equilibrium, with magnetic and thermal-pressure forces 
balancing the gravitational and external-pressure forces, in the absence of ambipolar diffu- 
sion. Any evolution at all is the result of the onset of ambipolar diffusion. The values of 
the free parameters are obtained from an initial reference state of the model cloud which is 
allowed to relax to an equilibrium configuration in the absence of ambipolar diffusion. The 
initial equilibrium state is thermally supercritical but magnetically subcritical. A detailed 
description of the calculation of the initial equilibrium state from a reference state is given 
in Ciolek & Mouschovias (1993). Magnetically critical or supercritical models have been 
studied by Fiedler & Mouschovias (1992) and Basu & Mouschovias (1995a, 1995b). 

At the outer boundary of the disk we require continuity of thermal as well as magnetic 
pressures 

p n (R,t) = ^, (20a) 
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and 

B z (R,t) = B ZjfSXt (R,t). 
On the axis of symmetry, radial components of all velocities vanish 
Vn (r = 0) = 0, v c {r = 0) = 0, v^r = 0) = 0, 



(20b) 



20c,d,e) 



v g 4r = 0) = 0, Vg+ (r = 0) = 0, v g0 (r = 0) = 0. ( 20f,g,h) 

In addition, the radial components of all forces vanish on the axis, which implies that 



and 



g r (r = 0) = 0, B r , z (r = 0) = 0, 



da n 



dr 



0, 



dB 



2,eq 



r=0 



dr 



0, 



r=0 



dz_ 

dr 



(20i,j) 



(20k,l) 



r=0 



6. THE CENTRAL-SINK APPROXIMATION 



In order to follow the evolution up to the time that 1 M© accumulates in a central region 
only a few AU in radius, Ri nner , some of the thin-disk equations must be modified and new 
equations must be added. The rate at which mass accumulates in this central region, is given 
by 

-27rr(7 n v n | r= „ (21a) 



dt 

Similarly, the rate of accretion of grain mass is 



g,* 



dt 

and that of magnetic flux 



= — 2nr (a. 



s- v g- 



--^innc 



9$ 



B,* 



dt 



-2vrr B z>eq v e \ R , 



(21b) 



(21c) 



where M*(t), M giJf (t) and $b,*(£) are the mass of gas, the mass of grains, and the magnetic 
flux in the central sink, respectively, destined to become a star. These three equations must 
be solved simultaneously with the other thin-disk equations. 

Also, the expressions for the radial component of the gravitational and magnetic fields 
must be modified in order to incorporate the influence of the mass and magnetic flux of the 
central cell, which, although not a part of the active computational region, affects the matter 
at r > -Ri nncr gravitationally and magnetically. No information about the spatial variation 
of the column density inside the central sink is available. However, we may account for the 
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formation of a dense core (at the center) at later times by expressing the gravitational field 



as 



= _uiiw + 2iiG / drV(7n (r')M(r,r'). (22) 
r Jo 

The first term on the RHS is the contribution of the mass in the central sink assuming that 
it is spherical in shape. In reality, the lower limit of the integral in the second term on the 
RHS of equation (22) should be -Ri nncr . However, for computational convenience we extend 
the integral to r = 0, and set a n (r < -Ri nncr ) = o"o> where a <C cr n (r = Ri miCT ), so as to avoid 
counting the mass inside the central sink twice. 

The radial component of the magnetic field on the surface of the disk, which appears in 
the magnetic tension force, is still obtained from the expression 



B 



poo 

(r) = - / dr'r'[B zm {r') ~ Bj[M{r, r'\ (23) 
Jo 



but -B 2jC q( r < -Rinncr) is now calculated by assuming that the flux inside the central sink 
is uniformly distributed. (Ciolek & Konigl 1998 used eq. [23] with the lower limit of the 
integral being -Rmner and they added a monopole term on the right-hand side to mimic the 
effect of the magnetic flux in the central sink.) 

The presence of a central point-like mass affects the equilibrium along the magnetic flux 
tubes as well. The equation of force balance in the z-direction now becomes 



,2 _ r , . K 2/ ,x . GM cent p n 



p Q (r,t)C 2 = P cxt + -Gaf 1 (r,t) + 



2 il \ i / 
r 



1 + 



<7 r 



2p Q r 



-1/2' 



(24) 



The last term on the right-hand side is the tidal gravitational stress corresponding to the z- 
component of the gravitational field of the central point mass (Ciolek & Konigl 1998). Near 
the central sink, tidal squeezing of the disk results in Z/r = a n /2p n < 1. So we may expand 
the last term and find that equation (24) becomes 



p n (r, t)C 2 = P ext + ^Gal(r, t) + (25 ) 



In order to follow shocks that might arise during the evolution, tensor artificial viscosity 
(Tscharnuter & Winkler 1979) is introduced. The viscosity tensor Q is given by 

5 = l 2 Pll (V-v)[Vv-^(V-v)l], tfV-v<0, 

= 0, if V • v > 0, (26) 
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where / is a length comparable to the local mesh spacing, and / is the unit tensor. In regions 
where the divergence of the velocity is negative, the radial viscosity force (per unit volume) 
in axisymmetric geometry is given by (Fiedler & Mouschovias 1992) 



fQ,r = -Ql 



r or 



dz r ' 



which, integrated over z in the thin-disk approximation, becomes 

2 1 o2 

1 On^n,r 0*V n 

3 



fQ,r = I 2 



1 v n , r da n dv Qtr 4 <7 n ( dv n ^ r 
3 r 



3 r dr dr 
4 <9t> n>r d 2 v ntr 



2 (9o- n / dv n r 



dr j > / 

2 



"3 an 



dr 2 
dv n . 



3 <9r 



dr dr 2 3 r 3 n,r 3 r 2 (9r r 2 n ' r (9r 



9r 



(27) 



(28) 



This term is added to the RHS of the radial force equation (18e). 



6.1. Initial and Boundary Conditions 

The calculation with the use of a central sink is considered a continuation of the previous 
simulation of the formation and contraction of a protostellar fragment in the following sense. 
When the central density exceeds ~ 10 11 cm" 3 and the assumption of isothermality is no 
longer valid, the simulation is stoped and the central region of radius i?i nn er is removed from 
the active computational mesh; it is now considered a sink cell. The calculation is then 
restarted, but the evolution is followed in detail only in the region r > -Ri nncr , in which 
the density is smaller than 10 11 cm" 3 . Thus the physical quantities at the end of the first 
calculation are used as initial data for the new calculation. 

At the inner boundary i?i nn er we impose perfectly absorbing boundary conditions by 
setting the values of all quantities at the boundary equal to those at the first mass zone 
(beyond the boundary) 

v n (r = 0) = v n (l), v e (r = 0) = v e (l), Vi (r = 0) = Vi (l), (29a,b,c) 

v s 4r = 0) = v g 4l), v g+ (r = 0)=v g+ (l), v g0 (r = 0) = Ug0 (l), (29d,e,f) 

<7 n (r = 0)=<7 n (l), <r g _(r = 0) = <7 g _(l), (30a,b) 

°g+( r = 0) = <7 g +(l), a g0 (r = 0) = <7 g0 (l). (30c,d) 

(Ciolek & Konigl 1998 used linear extrapolation of these quantities from the first computa- 
tional cell to the inner boundary.) The spatial derivatives at the first computational zone 
are calculated using one-sided differences instead of the central differencing used in the rest 



of the computational grid. The outer-boundary conditions (at r = R) are as in equations 
(20a) and (20b). 



7. THE DIMENSIONLESS PROBLEM AND METHOD OF SOLUTION 

We put the equations in dimensionless form by choosing units natural to the model 
molecular clouds. The unit of velocity is the isothermal sound speed in the neutral fluid, C; 
the units of column density <7 c ,rcf and acceleration 2nGa c ^ are those of the neutral column 
density and gravitational acceleration on the axis of symmetry of the reference state from 
which the initial equilibrium state is calculated. The unit of magnetic field is its uniform 
value in the reference state, B ref . The units of all other quantities are derivable from these 
four. 

We denote the dimensionless form of the cylindrical polar coordinates (r, z) by (£, () 
and the dimensionless time by r, so that equations (18a)-(18j) become 



<7 n (£,r) = 2p n (£,T)Z(£,T), 

d°n 1 <9(£q n ^n,g) n 
dr + £ d£ 

d(x g a Q ) 1 d[£<T n (xg-v g -A + Xg+v g +£ + Xgo%u)] 



dr 



Pn(£,r) 



1 



dr 



t d£ ctt d£ 

(3P cx t + el 



= 0, 



fa d(Tn _l j i? 



Pext + ^n(e,r)| , 



C cS = a 



mag,g 2 

Pd,cO 



C,eq 



C,eq 



+ 29e ; 


oo 

dee 




2 d 


' 1 




.i> 



3 2 



dZ 



poo 

buo = - / dee[B ( , cq (e) - b^m^ a 

Jo 



(31a) 
(31b) 

(31c) 

(31d) 

(31e) 

(31f) 



Bt.z^r ) + ( B^z-^r 



dZ 



-K 



(31g) 
(31h) 

(31i) 

(31j) 
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F T 

A mag ' ne 

c n A e 



0, 



V n + 



(31k) 

(311) 



Equations (181-t) are already in dimensionless form, provided that the velocities used are in 
units of C. After the introduction of the central sink, equations (21a)-(22), (24) and (28) 
become in dimenssionless form: 



dr 



9M g , 
dr 



^ ^inn 



dr 



s Jo 



Pn(£,r) = i(P cxt + ^)<;i + 



1/2' 



(32a) 
(32b) 
(32c) 
(32d) 
(32e) 



1 v n da n dv n 



4 <7 n / <9l> n \ 2 1 CT n f n d 2 V n 2 <9<T n / <9t> n 



3 f 3 £ v 0f 



3 ( 9( 2 3 <9£ V <9£ 



4 <9t> n d 2 v„ 1 cr 



+ 



n 2 



(32f) 



Equations (31a-l), (181-t) and (32a-f) comprise the dimensionless system to be solved nu- 
merically. The dependent variables are now dimensionless but in this section only, for ease 
in identification, we have retained the same symbols as for their dimensional counterparts. 
The quantities @ s and in the expressions for the velocities are already dimensionless 
(see Appendix B). The free parameters of the problem are: the initial central mass-to-flux 
ratio /id,co — (27rGr 1 / 2 cr Cjre f)/.B r . e / in units of its critical value for gravitational collapse in 
disk geometry, 1/(2ttG 1 ^ 2 ) (Nakano & Nakamura 1978); and the constant external pressure 
relative to the central gravitational "pressure" along the axis of symmetry of the reference 
state P C xt = Pext/(jG<j Ctie t). The reaction rates of the chemical network and the value of 
the radius of the dust grains are given in Appendix A. [The thermal (Jeans) length does not 
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appear explicitly in the equations because of the choice of the units we are using: The unit 
of length C 2 /27rG<7 C)re f is proportional to the Jeans length.] 

A control- volume method is employed for solving the dimensionless equations that gov- 
ern the evolution of the model cloud. The method consists of an implicit time integrator; 
a conservative upwind advective difference scheme that incorporates the second-order accu- 
rate monotonicity algorithm of Van Leer (1979); a second-order difference approximation of 
forces inside a mass zone; an integral approximation of the gravitational and magnetic fields; 
and an adaptive mesh capable of resolving the smallest natural lengthscale (i.e., thermal 
lengthscale At). The algebraic equations for the equilibrium abundances of charged species 
are solved at each time step. A detailed description of the numerical methods is given by 
Morton, Mouschovias k, Ciolek (1994). Once a central sink is introduced, we keep the grid 
stationary but nonuniform. 

8. DISCUSSION 

We have formulated the problem of the late accretion phase of a magnetic, isothermal 
disk surrounding a forming protostar. Our formulation is such that we can follow the struc- 
ture and evolution of the disk after the formation of an opaque core at the center of the disk, 
but before the protostar turns on and radiation from it affects the surrounding matter. The 
central-sink approximation allows us to follow the evolution of the system at significantly 
later times without the need of a numerically costly radiative transfer calculation for the 
opaque core. 

This is a similar approach to that of Ciolek and Konigl (1998) with two major physical 
and one numerical differences: (i) We have accounted for an additional charged species 
(positively-charged grains, which become important above densities ~ 10 9 cm -3 ), (ii) We 
have assumed that the magnetic flux is frozen only in the electron fluid and we have allowed 
the physics of the problem to determine the degree of attachment of the ions to the magnetic 
field lines. Detachment of ions in fact begins at densities ~ 10 8 cm -3 , as shown by Desch & 
Mouschovias (2001). (iii) We have used a more accurate calculation of the magnetic field and 
more appropriate boundary conditions. We have also included artificial viscosity in our code 
in order to resolve and track the formation and propagation of shocks. These differences 
result in new physical phenomena, such as the quasi-periodic formation and dissipation of 
shocks. The results of our numerical calculation are described in Paper II. 
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A. 



ABUNDANCES OF RELEVANT SPECIES 



The abundances of all relevant species (electrons, molecular ions, atomic ions, charged 
and neutral grains) are determined by assuming that chemical equilibrium is established and 
maintained throughout the evolution of the model molecular clouds. This assumption is 
valid since the chemical-reaction timescales are smaller than the evolutionary timescales of 
the model cores for neutral densities up to 10 14 cm -3 . 

In equilibrium, the rate of production of a specific species is equal to the rate of its 
destruction through the relevant chemical reactions listed in Table 1. Then the relative 
abundances of all species can be found solving the system of the creation/destruction rate 
equations for electrons, molecular ions, atomic ions, and positive grains: 



Since the total relative abundance of grains Xg an d the neutral density n n are obtained 
at each timestep from the continuity equations (31b) and (31c), the abundances of both 
negative and neutral grains can be determined from the system (Al)-(A5) once local charge 
neutrality is imposed: 



The relative abundance of species s is denoted by Xs(= n s/ n n)- 

The values of the rate constants for the gas phase reactions are taken from Ciolek & 
Mouschovias (1995) and references therein. For the dissociative recombination of HCO + , 
ctdr — 10~ 6 cm 3 s _1 ; for the radiative recombination of atomic ions, a TT = 10~ n cm 3 s _1 ; and 
for the charge exchange between atomic and molecular ions, /3 — 1.1 x 10~ 9 cm 3 s -1 . The 
typical value of the cosmic-ray ionization rate is (cr = 5 x 1CT 17 s -1 . 

The rate constants for electron-grain and ion-grain collisions (in cm 3 s _1 ) were originally 
calculated by Spitzer (1941, 1947) and later modified by Draine and Sutin (1987) to include 



CcR^n 
CcR^n 



(«drXm + + "rrXA+ + "egOXg + «eg+ Xg+ ) Xe™n > ( A l) 
("drXc + PXA° + "m+g-Xg- + "m+gOXgOXm+^n' ( A2 ) 

[«rrXc + "A+g-Xg- + a A +g°Xg°}XA+nl, (A3) 
(tteg+Xe + a g+g -x s -)Xg+nl, (A4) 
Xg-+Xg++X g o- (A5) 



= Xm+ + XA+ + Xg+ - Xe - Xg- ■ 



(A6) 
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corrections due to polarization of a grain by the approaching charged particle: 



Q; eg o = na 




(A7) 



na 



8/crT 



7rmi 



1 + 



7re^ 



2ak B T 



(A8) 



7ra 



nm 



ak B T 



1 + 



2 + (ak B T/e 2 ) 



(A9) 



= na 



8k B T 
nrrii 



1 + 



ak B T 



1 + 



2 + (ak B T/e 2 ) 



(A10) 



The quantities a eg o, a eg + are the rates of capture of e~ onto neutral and positive grains, 
respectively; and a ig o and a ig - are the rates of attachment of ions onto neutral and negative 
grains, respectively. Similarly, the rate of charge transfer and neutralization by direct grain- 
grain collisions is 



a 



g g 1 



n(a 1 + a 2 ) 



2 i 8k B T 



7im ved 



1 + 



ak B T 



1 + 



2 + (ak B T/e 2 



(All) 



In equations (A7)-(A11) e is the charge of the electron, m e and m\ the masses of the electron 
and ion, respectively, k B the Boltzmann constant, T the temperature, V e (V\) the sticking 
probability for an electron (ion) onto a grain, and m re d = miwi2/(mi + ^2) the reduced 
mass of the two grains. The terms enclosed in curly brackets account for the polarization of 
a grain by the incoming charge. 

The grains are assumed spherical with radius a = 3.75 x 10~ 6 cm and density p g = 2.3 g 
cm~ 3 , the average density of silicates. The temperature is kept at the constant value T = 10 
K. The mass of molecular ions is that of HCO, m m+ = 29 m p , while for atomic ions we use 
a mean value = 23.5 m p , in-between the mass of Na (m^ a = 23 m p ) and the mass of 
Mg (myig = 24 m p ). Since the masses of all the dominant ions are comparable, the fact that 
different ionic species dominate in different density regimes does not affect the evolution of 
the cloud cores. Finally, the sticking probability of electrons and ions onto grains are taken 
to be V c = 0.6 and V\ = 1, respectively (Draine & Sutin 1987; Umebayashi & Nakano 1980; 
Desch & Mouschovias 2001). 
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B. SPECIES VELOCITIES IN TERMS OF THE ELECTRON AND 

NEUTRAL VELOCITIES 

Once the velocity of the neutrals v n is known, the velocities of the other species can 
be found by solving the linear algebraic system of the force equations of the species and 
Ampere's law (eqs. [ld]-[lh] and [li], where j in eq. [li] is given by eq. [lj]). 

We add together equations (ld)-(lh) and we use equation (li) to eliminate j in the 
resulting equation to find the plasma force equation: 

Pn f \ . Pn / \ . Pn / \ . Pn / \ . Pn / v (V X B) X B 
(«e-«n)H («i-«n)H {v g --V D )-\ (f g+ - f n ) H {Vgp-V n ) = . 

T nc T n ; T~ng— 7~ng+ ^ngO 

(Bl) 

The terms on the left-hand side are the collisional coupling (momentum exchange) forces 
between the neutrals and the other five species (e, i, g_, g + and go). They were present 
on the right-hand side of the force equation (lc) of the neutrals, but were eliminated by 
using equation (Bl), which is a quantitative description of why and how the magnetic force 
appears in the force equation of the neutrals. 

We now use equation (15) for the electric field to eliminate E from equations (le)-(lg) 
and find that 

= + (B2) 

= _£!V( Vg __ Ve )xB + ^(« n - Vg _) + ^(« g0 -« g -) (B3) 

C " T gn T_ 

en g+/„. „ \ „ r» , Pg+/- - \ , PgO, 



= -^(v s+ -v e )xB + ^(v n -v g+ ) + ^(v s0 -v s+ ) (B4) 



Tgn 



where 



(^ + ^^ + ^^^) (B5) 

,7g0e,incl PgO T"g— i,inel PgO Tg-g+.incl, 

r = i^ + ^±^^ + ^±^^) {BG ) 



^TgOi.incl PgO Tg+e^nel PgO Tg+g-.inel, 

are the relevant timescales for inelastic collisions between grains. In particular, r_ is the 
timescale for a neutral grain to be involved in any inelastic reaction that converts negative 
grains to neutral grains or vice versa, and r + is the timescale for a neutral grain to participate 
in any inelastic reaction that transforms positive and neutral grains into each other. 

Equations (B1)-(B4), together with the neutral-grain force equation, 

= ^(v n - v s0 ) + ^(v g+ - v s0 ) + ^(« g _ - v s0 ) , (B7) 
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form a 10 x 10 linear system of equations, with unknowns being the r— and <p— components 
of the velocities of electrons, ions, and positive, negative, and neutral grains. 

We now integrate equations (B1)-(B4) and (B7) in the z-direction using the one-zone 
approximation (no variation of quantities with z inside the disk) and decompose them to r- 
and (^-components (assuming that there is no bulk rotation of the disk, v n ^ = 0). Note that 
ratios of mass densities, p, equal ratios of column densities, a, and the former rather than 
the latter are used whenever such a ratio appears. The resulting 10 x 10 linear system (with 
unknowns v h v- h<t> , v e , v e ^, v g -,v g - :l p, v g+ , v g+)<t> , v g0 , v g0 ,<t>) is: 

= WiTi n (^ - v ei(j >) + v n - Vi (B8) 
= uJiT in (v c - v{) - Vij (B9) 

= w g r gn (ue^-'u g -^)+u„-'Ug_ + -^- : ^(ugo-'Ug-) (B10) 

Pg- T - 

= LV g T gn (v g - - V e ) - U g _ ^ + ^(UgO,* - Ug-,*) (Bll) 

Pg- 

= uj g r gn (v g+>(j) -v c>(j> ) +v n -v g+ + -f^-^(v g0 -v g+ ) (B12) 

pg+ T+ 

= Lv g T gn (v e -v g+ )-v g+ ^ + ^^(v g0 , l p-v g+ , (t> ) (B13) 

pg+ T + 

= u n -Ugo + — (u g --Ugo) + — (ug+-Ugo) (B14) 

T_ T+ 
T T 

= -Ugo,* + — - Ugo,^) + — ( v g+,4> ~ v eP,4>) ( B15 ) 
= t^H ^i,«H v g +,*H fg-^H (B16) 

T n ; 7 ng+ ^ng— 7~ng0 

-fmagTtic . Tnc / \ . 7~ne / \ . 7~ne / \ . 7~ne / \ 

= V e -V n ^ (Wi-Wn)H K--^n)H (Wg + -W n )H (UgO-Un) 

<T n T n ; ^rig— 7 ng+ 7"ng0 

(B17) 

We first solve the system of equations (B8) - (B16) in that we express all the velocities in 
terms of v n and v e , and then substitute the resulting expressions for all velocities in equation 
(B17) to find v e as a function of v n . We note that in matrix form, the system of (B8) - (B16) 
is written as 

AV = v n C^ + v e C& , (B18) 
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where A is the 9x9 matrix of coefficients 



A 



-1 











^iTi„ 











-^iT in 


-WiTin 











1 

— 1 


n 
u 


n 
u 


n 
u 





n 
u 


— 1 — 1 — 


n 
u 


r_ 


n 
u 


-^ g T gn 


n 
u 


u 


CJgTg n 
















-1 -r_ 





r_ 











— 1 — r + 


r + 








a>gT gn 





— CUgTgn 








-cj g r gn 











— 1 — r + 


r+ 








T gn 


r gn 


T gn 




















t_ 



T+ 




TO 







T" gn 


Tg n 


_ r gn 

















^nc 


t_ 

Tne 


1 nc 


T 

Tnc 


1 


Tni 


T n g- 


T ng+ 


T ng 



with t = (1/Tgn + l/r- + l/r+) r_ = Pg T gn /(p g _T_), r + 
column vectors V, and are defined by 



(B19) 

Pgo r gn/(p g +T + ), and the three 



V 



V e ,<f, 



-1 


-1 


-1 


-1 









-^iT in 



CUgTg n 



— CJgTgj; 
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We use Cramer's method to solve the above system. We define 

D = det[A] . 



(B21) 



The notation D^ represents the determinant of A in which the column corresponding to the 
coefficients in front of the velocity of species s (where s is one of the following: i; i^; g_; g_^,; 
g + ; g + 0; go; go^; e^) has been substituted by C^ n \ Similarly, D* denotes the determinant of 
A in which the column s has been replaced by C^ e \ Then, the solution of the system is: 



Df Df 



D 



D 



v i,4> 



V g - = —^-Vn H 



D 



V n + 



D 



v g+ = 



D 
D n 
D 



-V n + 



D 
D 



D 



v g+,4> = 



D 



-Vn + 



Dl 



D 



D 



D 



(B22a,b) 
(B22c,d) 
(B22e,f) 



-26- 



S %n + Vg0> , = -g*Vn + -g**., (B22g,h) 



^ = + ^e- (B22i) 

Note that for all r-components of the velocities, D c s + = D. This can easily be seen if in 
the determinant D e s + £)" (which is the same as the determinant obtained by substituting 
the s-column by C e + C n ) we add to the s-column the 3 other columns that also correspond 
to r-velocities. The result is always the original determinant of the coefficients, D. This 
is a consequence of the property of the system of equations (B8) - (B17) that r-velocities, 
whether known or unknown, always appear as velocity differences in the equations. This is 
not the case for the ^-components of the velocities, where the assumption v n ^ = breaks 
the corresponding symmetry. 

This result is the mathematical expression of a physical property of the r-velocities, 
namely, that they vary smoothly between the r-velocity of the electrons, which are attached 
to the field lines, and the r-velocity of the neutrals. When a species is well attached to 
the field lines, its r-velocity coincides with that of the electrons. When collisions with the 
neutrals become frequent enough that the species is fully detached from the field lines, 
its r-velocity becomes equal to that of the neutrals. This physical picture can further be 
clarified by defining the quantities E s = {D* + D™)/D and 6 S = D e /D n , and noting that, as 
explained above, E s = 1 for all r-velocities. We can then rewrite the solution of the system 
as in equations (18a) - (18i). For the r-velocities S is then the attachment parameter (i.e. 
for Q s >■ 1, v s w v c , while for 9 S < 1, t) s « v n ). 

Finally, v e can be obtained as a function of v n by substituting equations (18a) - (18i) 
in equation (B17), which becomes 

-fmagTric / \ ( 1 7~nc @i Tnc ®g— ^~ne ®g+ ^~nc ®g0 



{ve _ Vn) M + + + + 



r ni 0; + 1 r ng _ g _ + 1 r ng+ 6 g+ + 1 r ng0 gO + 1 

(B23) 

Defining 

TV — 1 -\- Tnc ®' _|_ T nc Qg- _|_ r nc Qg+ _|_ r nc QgO (B24) 
r ni ©i + 1 Tng- ©g- + 1 r ng+ ©g+ + 1 T"ngO ©gO + 1 



we get equation (18k) which, if substituted in equations (18a) - (18i) gives the remaining 
r-velocities and the 0-velocities of all species as functions of v n . 
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Table 1: Chemical reaction network used in the calculation of the abundances of charged 
species. 

Relevant Chemical Reactions in Molecular Clouds 



Cosmic- Ray Ionization: H 2 + CR — > H^" + e 

H+ + H 2 -> H+ + H 

H+ + CO -> HCO+ + H 2 

Dissociative Recombination: HCO + + e — > H + CO 

Radiative recombination: A + + e — > A + /iz/ 

Charge tranfer: A + HCO+ -> A+ + HCO 

e _ attachment onto grains: e + go — > g_ 

e + g+ -»• go 

Atomic-ion attachment onto grains: A + + g_ — > A + g 

A + + go -> A + g + 

Molecular-ion attachment onto grains: HCO + + g — > HCO + g + 

HCO+ + g_ -> HCO + g 

Charge transfer by grain-grain collisions: g + + g_ — > 2g 



